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Abstract 

A  robust  method  for  computing  the  bed  shear  stress  in  unstratified 
combined  wave  and  current  flows  is  presented.  The  present  approach 
follows  from  existing  theories  describing  the  nonlinear  wave  and  current 
interaction  in  the  benthic  boundary  layer  but  is  designed  for  arbitrary 
wave,  current,  and  roughness  conditions,  including  the  limiting  case  of 
pure  waves  or  pure  currents.  The  stress  model  is  intended  as  a  stand-alone 
application  or  for  coupling  to  three-dimensional  shelf  circulation  models, 
where  a  broad  range  of  flow  conditions  are  encountered.  High-quality  data 
for  combined  flows  and  pure  waves  are  used  with  the  present  stress 
formulation  to  better  refine  empirical  model  closure  constants  in  the  fully 
rough  turbulent  regime.  Introducing  a  first-order  correction  to  the 
definition  of  the  wave  boundary  layer  thickness  produces  accurate 
estimates  of  both  the  measured  friction  factor  and  wave  boundary  layer 
height.  A  speed  of  convergence  test  indicates  that  the  present  model  is 
more  efficient  than  previous  models  that  use  the  same  turbulent  closure 
scheme.  This  is  primarily  due  to  an  improved  solution  algorithm  that 
avoids  the  nested  iterations  common  to  established  combined  wave  and 
current  bottom  boundary  layer  models. 
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1  Introduction 

An  important  physical  process  for  coastal  circulation  modeling  is  the 
interaction  in  the  bottom  boundary  layer  between  waves  and  currents  and 
how  these  both  interact  with  the  bottom  to  modify  bedforms  and  move 
sediment.  A  very  important  result  of  wave-current  interaction  theorized 
over  three  decades  ago  is  the  enhanced  current  shear  stress  due  to  waves 
(Smith  1975;  Grant  and  Madsen  1979).  Repeated  measurements  on  storm- 
dominated  shelves  have  illustrated  that  nonlinear  wave-current  interaction 
can  significantly  enhance  the  roughness  of  the  bed  and  the  stress  felt  by  the 
current  (e.g.,  Cacchione  and  Drake  1982;  Wiberg  and  Smith  1983;  Grant  et 
al.  1984;  Drake  et  al.  1992).  Therefore,  wave-current  interaction  is  expected 
to  play  a  dominant  role  in  the  momentum  balance  of  low-frequency  shelf 
motion  and  should  be  considered  in  any  realistic  modeling  effort  in  storm- 
dominated  shelf  regions.  This  is  especially  important  if  one  of  the  primary 
purposes  is  to  study  shallow  water  sediment  transport. 

1.1  Background 

Modeling  studies  of  shelf  circulation  patterns  that  incorporate  wave-current 
effects  in  the  bottom  boundary  layer  have  been  conducted  in  the  past  (e.g., 
Spaulding  and  Isaji  1985;  Cooper  and  Thompson  1989;  Signell  et  al.  1990; 
Davies  and  Lawrence  1994;  Keen  and  Slingerland  1993a,  b;  Keen  and  Glenn 
1994,  1995;  Warner  et  al.  2008;  Benetazzo  et  al.  2014).  Keen  and  Glenn 
(1994)  provide  a  brief  summary  of  coupled  and  uncoupled  versions  of  the 
Grant  and  Madsen  (1979)  (hereinafter  referred  to  as  GM)  and  Glenn  and 
Grant  (1987)  bottom  boundary  layer  models  (BBLMs)  implemented  in 
shelf-circulation  models.  Their  review  identifies  a  number  of  responses 
directly  related  to  enhanced  bottom  shear  stress  due  to  waves  on 
continental  shelves  including  a  reduction  in  current  speed  near  the  bottom, 
modification  of  sediment  transport  rates,  and  enhanced  turning  of  the 
current  vector  in  the  bottom  Ekman  layer  that  increases  coastal  upwelling 
and  downwelling.  Keen  and  Glenn  (1995)  also  showed  increased  offshore 
rotation  of  the  current  vector  during  downwelling  and  reduction  in  bottom 
current  speeds  in  shallow  water  in  a  simulation  of  storm  and  tidal  flow  in 
the  Middle  Atlantic  Bight.  Keen  and  Glenn  (1998)  carried  out  a  quantitative 
skill  assessment  of  model  performance  using  moored  current  meter  data 
from  the  Gulf  of  Mexico  during  Hurricane  Andrew.  One  of  the  sensitivities 
they  studied  included  a  three-order  of  magnitude  variation  in  bottom 
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roughness  length,  the  largest  roughness  serving  as  a  surrogate  for  the 
enhanced  apparent  bottom  roughness  known  to  occur  in  combined  wave 
and  current  flows.  Modeled  currents  showed  the  greatest  sensitivity  to  bed 
roughness  when  compared  to  bottom  currents  measured  in  a  water  depth  of 
15  meters  (m).  Normalized  peak  speed  differences  between  measured  and 
modeled  currents  decreased  when  the  apparent  roughness  was  increased 
from  0.1  centimeter  (cm)  to  10  cm.  Because  their  model  did  not  include 
wave-current  interaction,  the  roughness  and  stress  fields  could  not  evolve  in 
response  to  changing  wave  conditions.  Even  so,  the  higher  correlation 
between  modeled  and  measured  currents  in  shallow  water  for  simulated 
roughness  comparable  to  that  associated  with  the  presence  of  surface  waves 
reemphasizes  the  fact  that  wave-current  effects  are  important  on  storm- 
dominated  continental  shelves.  In  addition  to  these  earlier  studies,  there 
have  been  a  wide  variety  of  combined  flow  model  applications  in  shallow 
water  environments.  However,  development  of  new  analytical  models  has 
not  advanced  significantly  except  for  eddy  viscosity  formulations  (Yuan  and 
Madsen  2015;  Styles  and  Bryant  2016).  A  general  review  of  analytical 
combined  wave  and  current  models  is  provided  by  Nielsen  (1992)  and  Styles 
and  Bryant  (2016)  and  references  therein. 

The  above  results  are  based  on  a  streamlined  version  of  the  GM  wave  and 
current  BBLM  (Keen  and  Glenn  1994).  Like  the  original  GM  model,  the 
streamlined  version  assumes  that  the  roughness  length  is  small  compared 
to  the  wave  boundary  layer  height  and  that  the  height  of  the  reference 
current  needed  to  drive  the  model  (usually  the  lowest  grid  point)  is  greater 
than  the  wave  boundary  layer  height.  For  arbitrary  roughness  lengths  and 
model  grid  heights,  it  is  possible  that  under  some  conditions  neither  of 
these  requirements  will  be  met.  The  streamlined  version  also  uses  the 
discontinuous  eddy  viscosity  adopted  by  GM  and  Glenn  and  Grant  (1987), 
which  has  been  shown  to  be  less  accurate  than  more  physically  reasonable 
continuous  eddy  viscosity  profiles  for  combined  flows  (Glenn  1983;  Madsen 
and  Wikramanayake  1991;  Lynch  et  al.  1997;  Styles  and  Glenn  2000)  and 
pure  waves  (Sleath  1991;  Nielsen  1981, 1992;  Davies  and  Villaret  1997).  For 
this  study,  a  new  solution  to  the  Styles  and  Glenn  (2000)  BBLM  is  derived 
that  is  more  accurate  than  previous  models  referenced. 

1.2  Objective 

The  objective  of  this  technical  report  is  to  describe  the  theory  and 
equations  that  accompany  a  MATLAB  computer  program.  The  model  is  an 
extension  of  the  Styles  and  Glenn  (2000)  version  of  the  GM  model  but  has 
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been  modified  to  incorporate  arbitrary  roughness  configurations  and  a 
broader  range  of  turbulence  closure  schemes  (i.e.,  different  eddy  viscosity 
profiles).  A  secondary  objective  is  to  conduct  speed  of  convergence  tests  to 
gauge  model  performance  relative  to  other  widely-used  bottom  boundary 
layer  models. 

1.3  Approach 

The  approach  adopted  here  is  based  on  systematic  scaling  of  the  equations 
and  careful  selection  of  key  non-dimensional  parameters  to  produce  a 
closed  solution  that  is  numerically  stable.  The  algorithm  is  a  stand-alone 
application  that  can  be  used  to  predict  time-averaged  bed  shear  stress, 
maximum  bed  shear  stress  for  the  wave,  maximum  combined  bed  shear 
stress  for  the  wave  and  current,  bottom  roughness,  and  the  current  profile 
within  the  bottom  boundary  layer  for  combined  wave  and  current  flows.  It 
is  also  used  to  compute  the  bottom  stress  and  bed  roughness  in  three- 
dimensional  shelf  circulation  models  (Warner  et  al.  2008).  The  program 
also  predicts  sediment  transport  characteristics  for  non-cohesive  beds 
including  the  reference  concentration,  concentration  profile  within  the 
bottom  boundary  layer,  and  the  depth-integrated  concentration.  The 
sediment  transport  algorithms  can  further  be  used  to  predict  the 
concentration  profile  and  total  suspended  load  per  grid  cell. 

In  Chapter  2,  the  model  formulation  is  described  emphasizing  the 
modifications  required  to  extend  the  bottom  stress  theory  to  include  very 
rough  flow  conditions.  Model  sensitivities  to  the  eddy  viscosity  profile  and 
a  calibration  of  poorly  constrained  internal  model  closure  constants  is 
presented  in  Chapter  3.  This  is  followed  by  a  speed  of  convergence  test, 
and  the  results  are  summarized  in  Chapter  4. 
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2  Model  Formulation 


The  stress  model  developed  here  follows  that  of  GM,  in  which  the 
maximum  combined  shear  stress,  Tcw ,  is  written  as  the  vector  sum  of  the 
time-averaged  component  associated  with  the  current,  tc  ,  plus  the 
maximum  component  associated  with  the  wave,  Twm  , 


CW 


+ r 


wm’ 


(2-1) 


where  boldface  denotes  a  vector  quantity.  Writing  the  stresses  in  terms  of 
their  respective  shear  velocities,  u*  =  (r  /  p)1/2 ,  and  taking  the  magnitude 
gives 


/  2 
U 


2  u 


2  2 

lit 

c  wm 


COS  (j)Cy 


u. 


(2-2) 


where  p  is  the  fluid  density  and  </>cw(0  <tcw<ji/2)  is  the  angle  between 
the  wave  and  current.  To  obtain  a  closed  set  of  equations,  there  is  an 
adoption  of  the  usual  gradient  transport  relation  for  the  wave, 


u 


2 

*ivv 


du 

dz 


and  for  the  current, 


K 


dU 


dz 


(2-3) 


(2-4) 


where  K  is  the  time-independent  eddy  viscosity,  Uw  is  the  modulus  of  the 
wave  solution  in  the  lower  part  of  the  wave  boundary  layer,  U  is  the 
magnitude  of  the  horizontal  current,  z  is  the  vertical  coordinate  measured 
positive  upwards  from  the  bed,  and  z0  is  the  hydraulic  roughness.  Given 
profiles  for  the  eddy  viscosity,  wave,  and  current,  the  nonlinear  system 
Equations  (2-2),  (2-3),  and  (2-4)  can  be  solved  to  produce  the  bottom 
stress  vectors  Tcw ,  Twm ,  and  Tc . 

2.1  Small  to  intermediate  roughness 

For  conditions  in  which  the  height  of  the  roughness  elements  are  small  in 
comparison  to  the  wave  boundary  layer  thickness,  Glenn  (1983)  proposed 
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the  following  three-layer  continuous  eddy  viscosity  over  the  original  two- 
layer  discontinuous  formulation  used  by  GM: 

K  =  ku*cz  z>z2, 

K  =  ku*cwz  i  z1<z<z2,  (2-5) 

k  =  KILcwZ  Zq  <  z  <  z± , 


where  k  is  von  Karman’s  constant  (0.4),  z1  is  an  arbitrary  constant  scale 
height,  and  z2  =  z1  U*cw  /  U*c ,  which  is  determined  by  matching  the  eddy 
viscosities  at  z  =  z2  (Figure  1  [a]). 


Substituting  Equation  (2-5)  into  Equation  (2-3)  gives 


^ *wm  KU*cu)Zq  .  Tws, 

‘-cw 


where  the  non-dimensional  wave  shear  is  defined  by 


"P  —  lgv  ^U-iv 

ws~uu  dz 


1  du 


W 


uh 


f=?0 


(2-6) 


(2-7) 


and  ub  is  the  bottom  wave  orbital  velocity.  The  non-dimensional 
coordinate,  £  =  z  /  /CUI(£0  =  z0  /  lcw) ,  is  originally  derived  from  GM’s 
governing  equation  for  the  wave,  where  the  scale  height  of  the  wave 
boundary  layer  for  combined  flows,  lm ,  is  defined  by 

lcw=KU*cw/(0’  (2'8) 


and  0)  is  the  wave  radian  frequency.  Since  several  eddy  viscosities  will  be 
explored  in  this  analysis,  the  non-dimensional  wave  shear  is  introduced  as 
a  convenience.  The  solution  to  is  provided  in  Styles  and  Glenn  (2000). 
By  virtue  of  the  eddy  viscosity,  Equation  (2-6)  provides  a  relationship 
between  U*wm  and  U,cm . 
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Figure  1.  Schematic  illustrating  eddy  viscosity  profiles,  (a)  The  roughness  is  less  than  zi 
producing  the  three-layer  continuous  profile;  (b)  the  roughness  is  greater  than  zi  but  less  than 
zi,  so  the  eddy  viscosity  is  constant  near  the  bed;  (c)  the  roughness  length  is  large  compared 
to  the  wave  boundary  layer  height  and  increases  linearly,  indicative  of  a  pure  current 
boundary  layer;  (d)  the  eddy  viscosity  for  a  pure  wave  with  a  roughness  less  than  zi;  (e)  the 
roughness  is  greater  than  zi  for  a  pure  wave  producing  a  constant  eddy  viscosity;  and  (f)  the 
original  GM  discontinuous  eddy  viscosity  profile. 


Substituting  Equation  (2-5)  into  Equation  (2-4)  and  integrating  gives  the 
mean  current  profile: 


U{z) 

U(z) 

U(z) 


=  ^ln 

K 


+U(z2) 


U*c  zi\u{zA 


KU 


u ; 


-In 


KU, 


cw 


z2<z, 
z 1  z  <i  z2, 
Zq  z  , 


(2-9) 


where  the  no-slip  condition  at  z0  and  the  matching  requirement  that  the 
velocity  is  continuous  at  zx  and  z2  have  been  imposed.  Although 
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Equation  (2-9)  is  an  explicit  solution  for  the  current  in  this  application,  a 
current,  Ur ,  is  specified  at  a  given  height  above  the  bottom,  zr ,  and  u  c  is 
calculated  as  an  inverse  problem.  This  produces  a  relationship  between  U*c 
and  U*cw . 


2.2  Large  roughness 

If  the  roughness  length,  kb(=  30z0) ,  and  the  wave  boundary  layer 
thickness  are  of  the  same  order  of  magnitude,  then  the  boundary  layer  is 
fully  rough  turbulent,  and  it  becomes  possible  for  zo  to  exceed  lew.  Under 
these  conditions,  z0  can  become  greater  than  zi,  and  the  no-slip  condition  is 
applied  in  the  range  zi  <  z0  <  z2.  The  eddy  viscosity  is  constant  (Figure  1  [b]) 
near  the  bed,  and  the  profile  throughout  the  constant  stress  layer  is  given  by 


K  =  ku*cz  z2<z, 

K  =  ku*cwz1  z±<z0<z  <z2. 


(2-10) 


The  corresponding  kinematic  maximum  wave  stress  becomes 


H*u>m  KU*cu)Z±  7  > 

* cw 


(2-11) 


where  it  is  understood  that  Tws  implicitly  reflects  the  change  in  the 
solution  for  the  wave  shear  due  to  the  fact  that  the  eddy  viscosity  where 
the  no-slip  condition  is  applied  is  now  constant  instead  of  linearly 
increasing  as  in  Equation(2-5).  The  solution  for  Tws  is  presented  in 
Appendix  A.  Similarly,  the  solution  for  the  current  becomes 


U{z) 

U{z ) 


=  ^ln 

K 


+U{z2) 


_  U*c 


(z-z0) 


KU 


cw 


z2<z, 

Z^  Zq  z  <1  z2 , 


(2-12) 


where  the  no-slip  condition  is  now  applied  above  z1 . 

2.3  Large  roughness  with  vanishingly  small  waves 

For  the  case  of  the  vanishingly  small  waves,  z2  — >  z1  and  z0  may  even 
become  greater  than  z2 .  The  eddy  viscosity  where  the  no-slip  condition  is 
applied  then  becomes 
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K  =  KU*CZ  Z2<ZQ<Z  (2-13) 

(Figure  l  [c]),  so  that 

U*wm  KU,cZq  -[—Tws.  (2-14) 

Note  the  change  in  velocity  scale  from  U,.(W  to  U  (: ,  which  is  due  to  the  eddy 
viscosity  above  z2  being  a  function  of  U*c  and  not  U*cu; .  This  means  that  for 
very  small  though  finite  waves,  the  much  greater  shear  stress  associated 
with  the  current  can  still  affect  the  wave.  The  solution  using  the  eddy 
viscosity  given  by  Equation  (2-13)  is  in  Appendix  A.  Substituting  Equation 
(2-13)  into  Equation  (2-4),  the  current  reduces  to  the  classic  logarithmic 
profile, 


/  \ 


z2  Zq  <  z. 


(2-15) 


2.4  Non-dimensionalization  of  the  equations 

An  efficient  solution  algorithm  is  obtained  by  recasting  the  bottom  stress 
equations  into  a  suitable  non-dimensional  form.  Introducing  Wcw  as  the 
velocity  scale  produces  the  following  non-dimensional  parameters  that  will 
be  useful  in  formulating  the  bottom  stress  solution: 


cr 


_  ub 


_ 


U 


cw 


It 


II  —  ^ *wm 
/ ■*  . .  ’ 


cw 


U 


cw 


(2-16) 


where  a  is  related  to  a  combined  wave  and  current  friction  factor 
(a  =  l/4Zj2),  e  is  a  measure  of  the  relative  contribution  from  the 
current  to  the  total  stress,  and  ju  is  a  measure  of  the  relative  contribution 
from  the  wave  to  the  total  stress.  Squaring  both  sides  of  Equation  (2-2) 
and  rearranging  gives 

K  +  2u*wm  COS4X  +  Utwm  ~  ldcw  =  0,  (2-17) 

which  is  quadratic  in  u2c  with  the  solution 

K  =  ~u*wm  cos(/>cw  +  Vi4,m(cos2  (/lw  - 1)  +  idcw ,  (2-18) 
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where  the  positive  addition  of  the  second  term  is  chosen  to  ensure  that  u,c 
is  positive.  Dividing  both  sides  of  Equation  (2-18)  by  idcw  and  substituting 
ju  and  e  from  Equation  (2-16)  yields 


e2  =  -  n 2  cos(f>cw  +  ^4(cos2^cu;-1)  +  1, 
=  -H2  cos  (f)cw  +  ^/l  — jU4sin2  (fcw. 


The  kinematic  stress  for  the  wave  has  three  different  formulations 
corresponding  to  the  expressions  given  by  Equations  (2-6),  (2-11),  and 
(2-14).  Dividing  both  sides  of  these  equations  by  u2cw  and  substituting  // 
and  e  from  Equation  (2-16)  gives 

A*  = 

ILL2  =  K^oTws  £  < fo  < f2,  (2-20) 

l±2  =  K%2oYws  S2<  No¬ 
where  %1  =  z1/  lcw  and  £2  =  z2  /  Zcu) .  Similarly,  Equations  (2-9),  (2-12),  and 
(2-15)  can  be  used  to  formulate  three  separate  expressions  for  a .  Rather 
than  outlining  the  details  for  all  three  cases,  the  non-dimensionalization  is 
illustrated  using  Equation  (2-9).  Given  a  specified  current,  ur ,  at  a  height, 
zr ,  above  z2 ,  and  solving  for  o  yields 


(J 


KIL,. 


In 


“I- 1  —  e  +eln 


(2-21) 


If  the  observed  current,  ur ,  is  specified  at  values  of  zr  that  are  less  than  z2 
or  z1 ,  but  greater  than  z0 ,  three  more  equations  emerge.  The  resulting 
solutions  for  all  six  formulations  are  listed  in  Table  1,  along  with  their 
appropriate  ranges  of  validity. 

The  results  of  the  above  derivations  reveal  that  n  and  o  are  dependent  on 
non-dimensional  length  scales  that  arise  from  the  eddy  viscosity 
formulation  and  boundary  conditions.  These  unspecified  parameters  still 
must  be  determined  to  obtain  a  closed  solution.  In  general,  /i  is  a  function 
of  £0 ,  ,  <7 ,  and  e .  Only  the  first  two  variables  remain  unspecified,  and 

they  will  be  addressed  next. 
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Table  1.  Expressions  for  (7  derived  from  the  current  solution  discussed  in  the  text.  Inequalities  signify 


The  non-dimensional  roughness  height,  £0 ,  can  be  written 


where,  analogous  to  planetary  boundary  layers  (e.g.,  Grant  and  Madsen 
1986;  Wiberg  1995),  R  =  ilcw  /  z0co  is  an  internal  friction  Rossby  number 
for  combined  flows.  The  parameter  R  can  be  interpreted  as  the  ratio  of  the 
nonlinear  interaction  height  to  the  flow  roughness.  According  to  Madsen 
and  Wikramanayake  (1991),  the  dimensional  height  z1  is  a  function  of  the 
wave  boundary  layer  thickness,  z1  =  alcw ,  where  a  is  a  free  parameter  that 
represents  the  fraction  of  the  wave  boundary  layer  in  which  the  eddy 
viscosity  varies  linearly  with  height  and  is  determined  experimentally.  This 
defines  ^  =  a ,  which  gives  £2  =  a  /  e .  The  parameter  ju  is  now  a  function  of 
R,,  a,  a,  and  e . 

Examination  of  the  various  solutions  presented  in  Table  1  shows  that  at  a 
minimum 


0  =  0 


Uh  Z1  Zr 

— ,e,a 
Ur  Zo  Z\ 


(2-23) 


where  the  explicit  functional  dependence  on  z2  /  z1  and  zr  /  z2  has  been 
replaced  by  z2  =  z1  /  e .  Using  the  definition  for  lcw ,  z1  /  z0  can  be  written 
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z1  /z0  =  axR, .  An  analogous  expression  can  be  defined  for  z1/zr  (i.e., 
z1  /  zr  —  axR »  ),  where  R,r  =  u,CU)  /  (z^w) .  The  two  expressions  are  related 
by  R *  /  R*r  =  zr  /  z0 ,  where  zr  /  z0  is  an  independent  external  parameter. 

R*  and  cr  are  also  related  since  R,a  =  ub/  (oz0  =  Ab/z0,  where  Ab  is  the 
bottom  wave  excursion  amplitude.  Equation  (2-23)  is  now  an  implicit 
function  of  the  external  parameters  ub  /  ur,zr  /  z0,Ah  /  z0  and  the  internal 
parameters  a  and  e .  It  can  be  shown  that  e  is  a  function  of  Ab  /  z0,  a,  a , 
and  </>cw  so  that  cr  is  a  function  of  the  external  parameters 
ub/ur,  zr  /  z0,  Ab  /  z0,  <fcw ,  and  the  internal  closure  constant  a .  Although 
the  system  of  coupled  equations  resulting  from  this  analysis  does  not 
produce  an  algebraic  expression  for  the  shear  stresses,  a  closed  theoretical 
solution  exists.  The  nonlinear  system  therefore  can  be  solved  iteratively. 

2.5  Solution  algorithm  for  the  three-layer  model 

The  procedure  adopted  here  is  to  recast  the  series  of  non-dimensional 
expressions  derived  above  into  a  root-finding  algorithm  for  cr .  Applying  the 
pure  current  (u,  =  0)  or  pure  wave  (u  =  u  )  limit  shows  that  cr  is 
bounded  by  0  <  a  <  ub  /  u*wm .  The  solution  for  the  pure  wave  limit  is 
obtained  by  setting  /n  =  1  and  can  be  computed  independently  of  the 
combined  stress  solution.  Depending  on  the  root-finding  algorithm,  at  least 
one  initial  guess  for  cr  is  needed  to  start  the  iteration.  For  combined  flows, 
the  bisection  method  is  used  (Atkinson  1989)  since  the  root  is  guaranteed  to 
lie  between  the  pure  wave  and  pure  current  limits.  The  next  step  is  to 
determine  // ,  which  has  a  functional  dependence  that  can  be  described  by 


Ah 

zo 


(2-30) 


The  first  two  parameters  are  given,  the  third  is  assigned  an  initial  value  that 
lies  between  the  universal  limits,  and  the  last  parameter  is  unknown. 
Recalling  that  e  is  related  to  p  and  the  external  parameter  ^through 
Equation  (2-19),  an  internally  consistent  value  can  be  computed  by  recasting 
the  coupled  Equations  (2-19)  and  (2-20)  into  a  root-finding  algorithm  for  e 
similar  to  that  used  to  determine  o .  Again,  the  bisection  method  is  chosen 
since  e  is  bounded  by  universal  limits  (0  <  e  <  1) .  Once  p  and  e  have  been 
computed,  z1  /  z0  is  determined  by  z1  /  z0  =  axR  =  axAb  /  z0  /  o ,  and  z1  /  zr 
is  determined  by  z1  /  zr  =  axR  r  =  axAb  /  zr  /  o .  The  parameters  z2  /  zr  and 
z2  /  z0  are  related  to  z1  /  zr  and  z1  /  z0  through  e  ,  which  is  now  known. 
Given  the  above  estimates  for  these  non-dimensional  length  scales,  along 
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with  ub  /  ur ,  which  is  an  external  parameter,  a  new  value  for  o  is  computed 
from  the  equations  listed  in  Table  l.  The  relative  difference  between  the  new 
and  old  value  of  sigma  is  checked  to  see  if  it  is  below  some  prescribed 
tolerance.  If  it  is  not,  then  the  process  is  repeated  until  o  converges  (the 
tolerance  threshold  is  reached). 

2.6  Simplification  for  // 

The  solution  procedure  described  above  reveals  that  a  nested  iteration 
scheme  is  required,  in  which  an  inner  loop  is  first  initiated  to  produce 
internally  consistent  estimates  of  /ll  and  e  ,  and  then  an  outer  loop  is 
executed  to  solve  for  a .  These  iterations  represent  the  most 
computationally  expensive  operations  in  the  stress  solution.  If  the  inner 
loop  can  be  removed  from  the  solution  procedure,  then  the  total  number 
of  computations  will  be  reduced,  increasing  the  speed  of  convergence. 

An  examination  of  the  governing  equation  for  the  wave  (Styles  and  Glenn 
2000)  indicates  that  the  velocity  scale  (u*c)  for  the  stress  term  when  z  >  z2 
is  identical  to  the  formulation  above  the  wave  boundary  layer  derived  by 
GM.  Using  scaling  arguments  for  the  governing  equation  for  the  wave,  GM 
demonstrated  that  as  long  as  u»c  was  on  the  order  of  the  wave  velocity  or 
less,  then  the  stress  term  for  the  wave  outside  the  wave  boundary  layer 
could  be  neglected.  For  the  case  here,  which  considers  pure  currents  as  a 
possible  limit,  their  assumption  may  not  apply  when  the  current  is  much 
stronger  than  the  wave.  Under  these  circumstances,  the  wave  shear  and 
associated  wave  stress  for  z  >  z2  are  relatively  weak  so  that  the  wave 
solution  in  the  outer  region  is  well  described  by  the  linear  theory,  except 
possibly  under  very  rough  conditions  (Styles  and  Glenn  2000).  If  the 
stress  term  for  the  wave  is  neglected  above  z2 ,  then  the  solution  for  p 
becomes  independent  of  e  and  therefore  z2 .  This  eliminates  the  inner 
iteration  loop  required  to  produce  an  internally  consistent  value  for  ^  and 
e  and  accelerates  the  speed  of  convergence  without  appreciably  altering 
the  results  of  the  stress  model  based  on  a  three-layer  eddy  viscosity  for  n  . 

Appendix  B  presents  a  derivation  of  Tws  and  fi  based  on  a  simpler, 
continuous  two-layer  eddy  viscosity  (Figure  1  [d],  [e])  in  which  the  stress 
term  for  the  wave  above  z2  is  neglected.  The  solution  procedure  follows 
that  described  in  Section  2.5  except  that  fi  no  longer  depends  on  e  . 
Instead,  e  is  computed  explicitly  through  Equation  (2-19). 
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3  Model  Results 

To  illustrate  the  properties  of  the  stress  model,  results  are  presented  based 
on  the  simplified  solution  for  fi  discussed  in  Section  2.6.  The  input 
parameters  consist  of  the  external  variables  Ab  /  z0,  zr  /  z0,ub  /  ur ,  and  <f)cw 
and  the  internal  closure  constant  a . 

Because  the  stress  model  is  designed  for  a  broad  range  of  input  wave  and 
current  conditions  that  may  be  produced  by  a  shelf  circulation  model,  the 
parameter  ranges  for  Ab  /  z0  and  zr  /  z0  are  10  ~3<Ah/z0<  106  and 
1.01  <  z;,  /  z0  <  106 ,  respectively.  The  lower  limit  is  chosen  to  represent  a 
maximum  relative  roughness  for  the  wave  (z0  /  Ab)  of  103 .  The  upper 
limit  is  chosen  to  represent  a  too  cm  current  height  (zr)  or  excursion 
amplitude  ( Ab )  over  a  ripple-free  bed  with  a  minimum  grain  diameter 
roughness  of  ~io_3  cm  (10  microns  [pm]).  Other  model  parameters  have 
been  fixed  with  values  of  a  =  1,  <fcw  =  0 ,  and  uh/ur  =  1 .  The  latter  is 
chosen  so  that  the  wave  and  current  outside  the  wave  boundary  layer  are 
approximately  the  same  order  of  magnitude.  Past  expressions  for  a  have 
ranged  between  approximately  0.15  and  2  (Glenn  1983;  Madsen  and 
Wikramanayake  1991;  Lynch  et  al.  1997).  Therefore,  an  intermediate  value 
(1.0)  is  chosen  to  illustrate  the  model  characteristics. 

3.1  Fundamental  model  characteristics 

Figure  2  depicts  selected  internal  model  parameters  identified  in  the  text 
as  a  function  of  the  independent  external  parameters  zr  /  z0  and  Ah/  z0. 
Individual  model  parameters  show  varying  degrees  of  sensitivity  to  zr  /  z0 
and  Ab/  z0 ,  especially  for  extreme  values.  The  first  two  parameters, 

H  =  u,wm  /  u*cw  and  e  =  u,c  /  u,cw ,  illustrate  the  dynamic  features  of  the 
stress  model  since  they  define  the  relative  proportions  of  the  wave  and 
current  shear  velocities  to  the  total.  The  pure  wave  or  pure  current  limits 
are  easily  interpreted  graphically  as  /j  — *•  1  and  e  ->  0  for  pure  waves,  and 
ji  — ►  0  and  e  — >  1  for  pure  currents.  As  zr  /  z0  — >  1  for  a  fixed  ur ,  the 
current  shear  becomes  very  large.  This  produces  a  larger  stress  for  the 
current  (e  — >  1)  that  will  dominate  over  the  wave  in  the  combined  flow.  As 
this  ratio  increases,  the  current  shear  begins  to  decrease,  with  an 
associated  reduction  in  e  .  Eventually,  zr  /  z0  will  become  so  large  that  the 
turbulent  stresses  associated  with  the  current  (for  constant  ur )  must 
vanish  and  e  =  0 .  In  this  limit,  the  solution  becomes  that  of  a  pure  wave 
(/j  — s- 1) .  The  rate  at  which  the  pure  wave  limit  is  approached  is  also  a 
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function  of  Ah  /  z0 .  For  constant  ub ,  decreases  in  Ah  /  z0  lead  to  greater 
frictional  drag  for  the  wave  and  an  associated  increase  in  bottom  stress. 
This  behavior  is  apparent  in  the  first  two  plots,  as  the  pure  wave  limit 
proceeds  more  rapidly  as  a  function  of  zr  /  z0  for  smaller  Ab/  z0. 


Figure  2.  Selected  internal  model  parameters  as  a  function  of  zr  /  z0  and  Ab  /  z0 . 
Definitions  of  the  internal  variables  are  provided  in  the  text.  Values  for  Ab  /  z0  range 
from  lO3  to  106  in  decadal  increments. 
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Another  notable  feature  is  that  both  /./  and  e  become  independent  of 
Ab  /  z0  when  this  ratio  is  less  than  l.  This  can  be  understood  by  examining 
the  non-dimensional  length  scales  z1  /  z0  and  z2  /  z0 ,  and  the  parameter  o 
as  a  function  of  zr  /  z0 ,  when  10  3  <  A/z„<  10  1 .  Figure  2  reveals  that 
z1  /  z0  and  z2  /  z0  are  always  less  than  one  in  this  range  so  that  the  eddy 
viscosity  at  the  bed  becomes  constant  and  ju2  =  Tc^/qa  (Appendix  B). 
Substituting  o  from  Equation  (2-29)  into  the  above  expression  for  //  gives 


zr>z0>z2>z1,  (3-1) 


which  is  independent  of  Ab  /  z0 .  The  parameter  e  is  related  to  q  through 
Equation  (2-19)  and  is  also  independent  of  Ab/  z0 . 

The  second  group  of  parameters,  z1  /  zQ,z2  /  zQ,zx  /  zr,  and  z2  /  zr, 
represents  the  length  scales  of  the  flow.  All  four  parameters  exhibit  a 
strong  dependence  on  Ab/  z0,  but  only  z±  /  zr  and  z2  /  zr  are  dependent  of 
zr  /  z0  when  this  ratio  becomes  very  large.  An  important  consideration  for 
modeling  applications  is  that  all  four  parameters  are  well  behaved  for  the 
broad  range  of  Ab  /  z0  and  zr  /  z0  used  here.  The  only  exception  occurs 
when  zr  /  z0  =  1  (ur  =  0) ,  which  is  a  degenerate  case. 

Tws  is  sensitive  to  smaller  values  of  zr  /  z0  and  to  Ab  /  z0  as  long  as 
z1  /  z0  >  1 .  When  z1  /  z0  <  1 ,  the  eddy  viscosity  that  defines  the  maximum 
wave  shear  is  constant,  and  the  resulting  equation  for  Tws ,  which  is  derived 
in  Appendix  B,  is  independent  of  Ab/  z0 .  The  final  two  parameters,  R *  and 
o ,  are  both  sensitive  to  smaller  values  of  zr  /  z0 ,  and  R  is  sensitive  to 
Ab  /  z0  for  all  values.  The  parameter  o  is  not  very  sensitive  to  the  very  large 
changes  in  either  zr  /  z0  or  Ab  /  z0 .  This  demonstrates  the  advantage  of 
selecting  a  as  the  function  best  suited  for  a  root-finding  algorithm  to  close 
the  stress  solution. 


The  stress  model  also  depends  on  ub/ur,  which  is  a  measure  of  the  relative 
strength  of  the  wave  to  the  current.  Figure  3  shows  n,  e,  and  o  for  the 
same  conditions  illustrated  in  Figure  2  but  with  ub/ur=  10  (large  wave) 
and  0.1  (small  wave).  The  general  trends  are  the  same  as  the  ub  /  ur  =1 
case,  with  the  exception  that  for  larger  ub  /  ur  the  solution  approaches  that 
of  a  pure  wave  much  faster  as  a  function  of  zr  /  z0 ,  and  for  ub  /  ur  =  0.1  the 
change  is  more  gradual.  Since  a  larger  ub  /  ur  signifies  a  stronger  ambient 
wave,  it  is  expected  that  the  solution  should  converge  to  the  pure  wave 
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limit  much  more  rapidly  as  zr  /  z0  increases.  The  opposite  is  true  for  a 
relatively  large  current  ( ub  /  ur  =  o.i).  The  parameter  o  also  is  not  as 
sensitive  to  Ab  /  z0  when  the  current  is  much  stronger  than  the  wave. 


Figure  3.  Similar  to  Figure  2  but  showing  only  ju ,  e  ,  and  <7 .  Left  column  is  for 
Ub  /  Ur=  10  (large  waves),  and  right  column  is  for  Ub  /  Ur  =  0.1  (large 

currents). 


3.2  Sensitivity  to  the  direction  between  the  wave  and  current 

Equation  (2-19)  expresses  a  closed  relationship  between  e  and  p  given  the 
external  parameter  (f>cw .  A  plot  of  e  as  a  function  of  p  for  several  values  of 
(fcw  is  shown  in  Figure  4.  For  small  n ,  e  is  not  very  sensitive  to  n  or  (fc.w .  As 
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//  becomes  larger,  and  the  shear  stress  associated  with  the  wave  becomes 
significant,  e  becomes  much  more  sensitive  to  <f>cw .  This  sensitivity  can  be 
described  mathematically,  since  Equation  (2-19)  reduces  to 

e2  —1  —  /J2  (3-2) 

for  co-directional  flow  (</>cw  =  0)  and  to 

e2  =  Vl-U4  (3-3) 

for  orthogonal  flow  {(j)cw  —  n/2).  The  larger  exponent  for  n  in 
Equation  (3-3)  tends  to  make  e  larger  when  the  wave  and  current  are  at 
right  angles.  Physically,  this  means  that  for  a  fixed  maximum  wave  stress 
vector  in  the  presence  of  a  current,  the  magnitude  of  the  time-averaged 
shear  stress  must  continually  increase  as  (fcw  goes  from  o  to  n  /  2 ,  if  the 
magnitude  of  the  maximum  total  stress  vector  is  to  remain  constant.  The 
direction  of  the  maximum  total  stress  vector  will  of  course  change  as  the 
time-averaged  shear  stress  vector  rotates  toward  n/2.  During  storms,  the 
wave  and  current  vectors  near  the  coast  are  generally  at  a  high  angle 
(closer  to  orthogonal  than  co-directional)  and  both  are  relatively  strong.  If 
topographic  steering  or  an  evolving  current  (e.g.,  tides)  produces  local 
regions  where  <f)cw  becomes  small,  an  associated  increase  in  the  magnitude 
of  the  total  shear  stress  may  occur.  Several  studies  have  examined 
combined  flows  in  which  the  angle  between  the  wave  and  current  is  90 
degrees  (Fernando  et  al.  2011;  Madsen  et  al.  2008;  Musumeci  et  al.  2006). 
There  have  been  fewer  observational  studies  to  characterize  the  stress 
within  the  wave  boundary  layer  for  arbitrary  wave  and  current  vectors. 
Some  studies  indicate  that  a  first-order  effect  is  a  reduction  in  the  bottom 
roughness  for  the  current  as  <pcw  increases  (Sorenson  et  al.  1995;  Styles 
1998;  Madsen  et  al.  2008;  Styles  and  Bryant  2016). 
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Figure  4.  Sensitivity  of  the  parameters  e  and  //  to  the  direction 
between  the  wave  and  the  current. 


3.3  Bottom  stress  sensitivity  to  the  eddy  viscosity  profile 

The  solution  presented  above  neglects  the  stress  term  in  the  governing 
equation  for  the  wave  above  z2  but  retains  it  for  the  current.  This  was 
justified  based  on  the  assumption  that  u*c  was  on  the  order  of  the  wave 
velocity  or  less  (GM).  Styles  and  Glenn  (2000)  also  have  argued  that  the 
details  of  the  eddy  viscosity  outside  the  wave  boundary  layer  are  not 
important  in  determining  the  bed  stress  except  possibly  for  very  rough 
beds.  Both  the  Styles  and  Glenn  (2000)  three-layer  eddy  viscosity  depicted 
in  Figure  1  (a,  b,  c)  and  the  simplified,  two-layer  continuous  eddy  viscosity 
for  the  wave  depicted  in  Figure  1  (d,  e)  are  identical  as  z  — >•  z0  but  diverge 
above  z2 .  Another  eddy  viscosity  profile  that  has  been  used  extensively  in 
the  past  is  the  linearly  increasing,  discontinuous  form  originally  proposed 
by  GM  (Figure  1  [f]).  Since  all  three  formulations  are  different  above  the 
wave  boundary  layer  but  the  same  below  z1 ,  the  present  stress  model  can 
be  used  to  examine  how  the  details  of  the  eddy  viscosity  profile  outside  the 
wave  boundary  layer  affect  bed  stress  estimates. 

A  stress  model  based  on  the  GM  eddy  viscosity  does  not  include  the  z1  or  z2 
terms.  Instead,  GM  prescribe  the  height  of  the  wave  boundary  layer,  Scw , 
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which  is  also  formulated  as  a  constant,  n,  times  lcw  (8CW  =  nlcw).  For  this 
comparison,  n  is  set  equal  to  2,  which  is  the  typical  value  used  in 
applications  (Glenn  and  Grant  1987;  Madsen  et  al.  1993;  Madsen  1994; 
Keen  and  Glenn  1994).  To  highlight  the  differences  between  the  three  eddy 
viscosity  formulations,  the  ranges  of  the  input  variables  are  reduced  to 
10  1  <Ab/z0<  103  and  1.01  <zr  /  z0  <  103 .  Other  parameters  are  set  with 
values  of  0cw  =  O,ub/ur=l  and  a  =  1 .  To  illustrate  the  dynamical 
properties  of  the  three  modeling  approaches,  the  comparison  focuses  on  the 
parameters  p,  e,  and  a ,  which  are  depicted  in  Figure  5.  For  Ahj  z0>  100 , 
all  three  models  produce  approximately  the  same  result,  and  therefore,  the 
solution  is  not  sensitive  to  the  form  of  the  eddy  viscosity  in  the  outer  wave 
boundary  layer  and  above.  This  point  is  argued  by  Styles  and  Glenn  (2000), 
who  claim  that  the  stratification  correction  also  introduces  arbitrary 
changes  to  the  eddy  viscosity  and  therefore  can  be  neglected  in  the  wave 
stress  solution. 


Figure  5.  Comparison  of  /i,  e, ,  and  <7  for  the  GM  (thin  solid).  Styles  and  Glenn  (2000)  (dash),  and 
two-layer  continuous  (thick  solid)  eddy  viscosity  profiles.  Ab  )  ZQ  ranges  from  10 1  to  103  in  decadal 

increments. 


For  a  larger  relative  roughness,  the  three  solutions  begin  to  diverge.  This  is 
most  apparent  for  the  two-layer  continuous  eddy  viscosity,  which  was 
shown  to  produce  an  upper  bound  on  q  and  a  lower  bound  on  e  and  a 
when  z1  /  z0  <  1 .  The  GM  and  Styles  and  Glenn  (2000)  eddy  viscosity 
profiles  for  the  wave  are  not  constant  in  the  outer  portion  of  the  constant 
stress  layer  so  that  a,  ju,  and  e  remain  functions  of  Ab  / z0  when  the 
relative  roughness  is  very  large. 

In  the  limit  of  a  pure  current  (zr  /  z0  -» 1) ,  both  n  and  e  based  on  the 
stress  models  derived  using  the  GM  and  Styles  and  Glenn  (2000)  eddy 
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viscosities  converge.  This  limit  was  mentioned  above  as  a  possible  case 
when  uc  might  be  important  in  the  governing  equation  for  the  wave 
outside  the  wave  boundary  layer.  Since  the  result  based  on  the  Styles  and 
Glenn  (2000)  eddy  viscosity  profile  includes  the  stress  term  above  z2  and 
the  GM  eddy  viscosity  does  not,  the  importance  of  retaining  the  stress 
term  for  this  case  appears  minimal.  The  fact  that  the  three  models  produce 
dissimilar  results  for  a  large  relative  roughness  (Ab  /  z0  <  10)  suggests  that 
careful  consideration  of  the  parameterization  of  the  eddy  viscosity  for  a 
rough  bed  is  important. 

3.4  Evaluation  of  a 

The  above  analysis  has  revealed  that  the  model  is  most  sensitive  to  the  eddy 
viscosity  profile  for  very  rough  beds  ( kb  /  Ah  >  1.0) .  On  storm-dominated 
sandy  continental  shelves,  the  roughest  beds  are  speculated  to  be  associated 
with  the  presence  of  relic  ripples,  which  can  have  maximum  ripple  heights 
that  exceed  10  cm  (Traykovski  et  al.  1999).  The  amount  of  time  that  relic 
ripples  dominate  the  roughness  signature  on  sandy  continental  shelves  is  a 
function  of  the  state  of  ripple  degradation  after  storms  (Nelson  and 
Voulgaris  2014).  Assuming  that  relic  ripples  persist  for  some  time  after 
storm  events,  it  is  possible  to  estimate  the  average  amount  of  time  that  relic 
ripples  may  be  present.  Studies  of  storm-forced  transport  on  the  New  Jersey 
shelf  (Styles  1998)  have  indicated  that  the  average  storm,  as  defined  by  the 
time  that  the  shear  stress  based  on  skin  friction  exceeds  the  minimum  for 
the  initiation  of  sediment  motion,  lasts  approximately  24  hours  and  that 
approximately  10  such  storms  occur  annually.  Assuming  that  biological 
activity  sufficiently  degrades  ripples  within  a  week  or  two  after  a  storm 
(Traykovski  et  al.  1999)  gives  an  annual  maximum  relic  ripple  period  of  2  to 
5  months.  This  can  be  a  significant  amount  of  time  and  indicates  that 
BBLMs  must  be  designed  to  accommodate  a  variety  of  roughness  conditions 
that  are  likely  to  occur  on  wave-dominated  continental  shelves. 

Model  sensitivity  to  the  form  of  the  eddy  viscosity  profile  can  be 
investigated  by  modulating  the  parameter  a .  If  a  is  relatively  large,  then 
the  eddy  viscosity  in  the  vicinity  of  z0  increases  linearly  with  height.  This 
is  the  same  as  the  GM  model  deep  within  the  wave  boundary  layer.  For 
intermediate  values  of  a ,  the  eddy  viscosity  is  identical  to  Madsen  and 
Wikramanayake  (1991)  and  Styles  and  Glenn  (2000).  If  a  is  very  small, 
then  the  eddy  viscosity  in  the  vicinity  of  z0  is  constant,  and  it  is  similar  to 
vertically  uniform  profiles  that  have  been  suggested  for  very  rough  beds 
(Nielsen  1992;  Sleath  1991).  Modifying  a  in  the  present  model  effectively 
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reproduces  the  eddy  viscosity  profiles  discussed  above,  which  were  shown 
to  produce  different  results  when  kb  /  Ah  was  greater  than  approximately  l. 

For  a  given  set  of  external  wave,  current,  and  roughness  conditions,  the 
boundary  shear  stress  becomes  sensitive  only  to  the  value  of  a .  In 
laboratory  flumes,  all  of  these  parameters  can  easily  be  prescribed  or 
measured  independently  of  a  combined  flow  model  with  the  exception  of 
the  bottom  roughness.  This  is  because  the  precise  mathematical 
formulation  depends  on  the  size  and  shape  of  the  bedforms  present,  which 
means  that  kb  varies  as  a  function  of  the  experimental  conditions  and  is 
not  universal  in  form.  Therefore,  the  experimental  setting  must  conform 
as  closely  as  possible  to  the  actual  environmental  conditions  to  which  the 
calibration  results  apply.  In  this  case,  the  experimental  conditions  must 
include  roughness  elements  that  simulate  the  approximate  shape  of  wave¬ 
generated  ripples  and  more  importantly  they  must  return  a  consistent 
roughness  value  based  on  several  independent  methods  of  determination. 

The  data  sets  used  to  evaluate  a  for  combined  flows  are  obtained  from 
Mathisen  and  Madsen  (1996a,  b)  (hereinafter  referred  to  as  MM).  MM 
conducted  detailed  experiments  of  co-directional  wave  and  current  flows 
in  a  laboratory  flume.  To  ensure  rough  turbulent  conditions,  they  modified 
the  bed  of  their  flume  with  triangular  shaped  bars  that  were  scaled  to 
simulate  the  geometry  of  two-dimensional  (2-D)  wave-generated  ripples. 
In  nearly  all  of  their  experiments,  kb  /  Ab  >  1 ,  which  is  ideal  for  evaluating 
the  present  stress  model  for  rough  conditions.  MM  reported  all  necessary 
input  data  to  drive  the  model  including  bottom  roughness  height  and  the 
wave  friction  factor,  which  was  determined  by  measuring  the  decay  in 
wave  height  over  the  length  of  the  flume  and  relating  that  to  dissipation 
due  to  bottom  friction.  Noting  that  (1/  a  —  yjfw  /  2  ,  friction  factor  curves 
can  be  generated  and  compared  to  their  measurements.  Setting  kb  equal  to 
the  roughness  determined  for  pure  currents  (MM  1996a,  Table  2)  and 
using  the  measured  water  particle  amplitudes,  orbital  velocities  and  mean 
currents,  a  family  of  friction  factor  curves  as  a  function  of  a  is  generated. 
The  current  roughness  is  chosen  since  MM’s  results  demonstrated  that  kb 
was  nearly  the  same  for  waves  in  the  presence  and  absence  of  currents, 
and  currents  in  the  presence  and  absence  of  waves.  Also,  their  current 
roughness  estimates  are  independent  of  the  GM  combined  wave  and 
current  model,  whereas  this  model  is  used  to  determine  the  roughness  for 
all  their  cases  with  waves.  As  noted  above,  increasing  a  leads  to  an  eddy 
viscosity  profile  that  is  very  similar  to  the  GM  formulation.  The  decision  to 
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use  the  roughness  length  for  pure  currents  ensures  that  the  present 
method  to  determine  a  is  not  inherently  dependent  on  the  GM  model 
(through  the  bottom  roughness),  which  may  bias  the  results  to  favor  larger 
values  of  a .  To  quantify  the  comparison  between  the  model  and  data, 
there  is  an  adoption  of  the  relative  error  defined  by 


ln(e) 


1 

4Efln(y)-in(y))2f 

1V  f=i 


(3-4) 


where  Y{  is  the  measured  data  point,  F  is  the  corresponding  model 
estimate,  and  N  is  the  number  of  data  points  (Wikramanayake  and 
Madsen  1991).  The  friction  factor  curve  that  minimizes  e  identifies  the 
optimum  a . 

Using  values  that  range  from  0.15  to  2  (Glenn  1983;  Madsen  and 
Wikramanayake  1991;  Lynch  et  al.  1997;  Styles  1998),  the  lowest  error  (e  = 
1.3)  is  obtained  when  a  is  set  equal  to  0.75.  The  corresponding  modeled 
wave  boundary  layer  thickness,  as  determined  by  8CW  =  nlcw  (n  =  2) ,  is  a 
factor  of  3  too  low  when  compared  to  the  observed  height  derived  from 
wave  and  current  profile  measurements  in  the  flume.  MM  also  under¬ 
estimate  the  height  of  the  wave  boundary  layer  using  the  GM  model.  MM 
attribute  the  enhanced  boundary  layer  thickness  to  the  increased  bed 
roughness  associated  with  their  fixed  artificial  roughness  elements.  Relic 
ripples  may  play  a  role  similar  to  artificial  roughness  elements,  and 
therefore,  may  produce  an  enhanced  boundary  layer  thickness  relative  to 
smoother  flow  conditions.  Although  the  main  purpose  here  is  to  describe  an 
algorithm  to  compute  the  total  bed  stress,  the  resulting  wave  and  current 
stress  components  are  integral  components  of  suspended  sediment 
concentration  and  velocity  profile  models  (i.e.,  Smith  1975;  Wiberg  and 
Smith  1983;  GM;  Glenn  and  Grant  1987).  A  model  designed  to  predict 
current  and  suspended  sediment  concentration  profiles  should  be  able  to 
reproduce  accurately  both  the  friction  factor  (stress)  and  the  wave  boundary 
layer  thickness. 

The  fact  that  the  model  tends  to  underestimate  the  thickness  of  the  wave 
boundary  layer,  yet  accurately  predicts  the  wave  friction  factor,  suggests 
that  it  is  the  internal  length  scales,  which  define  the  height  and  thickness 
of  the  various  regions  in  the  boundary  layer,  as  opposed  to  the  velocity 
scales,  which  define  the  shear  stresses,  that  should  be  reexamined.  The 
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present  and  GM  model  have  adjustable  internal  length  scales  z1  and  8CW , 
respectively.  It  is  hypothesized  that  for  very  rough  conditions  ( kb  /  Ab  >  1) , 
the  constants  multiplying  lcw  in  the  definition  of  z1  and  8CW  are  now 
functions  of  the  relative  roughness.  A  very  simple  approximation  that 
incorporates  explicitly  the  relative  roughness  in  the  definition  of  z1  but 
reverts  to  the  existing  formulation  in  the  limit  kb  /  Ab  — >■  0  is  to  modify  z1 
as  z1  =  alcw(a0  +  ajch  /  Ab  +  a2(kh  /  Ah)2  +  Since  kb  /  Ab  is  expected  to 
become  a  leading  order  term  only  for  very  rough  beds,  the  series  can  be 
truncated  to  first  order  to  give  z1  =  crfcw(l  +  j 8kb  /  Ab) .  A  similar  expression 
is  proposed  for  8CW :  8CW  =  nlcw{  1  +  /3kb  /  Ah) .  For  smooth  to  moderately 
rough  turbulent  conditions  (kb  /  Ah  <  0.1) ,  the  GM  model  has  been  shown 
to  accurately  predict  the  shear  stress  and  apparent  roughness  with  n  =  2 
(e.g.,  Grant  et  al.  1984;  Drake  and  Cacchione  1992;  Drake  et  al.  1992).  This 
leaves  two  undetermined  parameters  ( a  and  (3 )  that  must  be  calibrated 
from  data.  Optimal  values  are  found  by  choosing  the  combination  that 
minimizes  the  relative  difference  for  both  the  friction  factor  and  the  wave 
boundary  layer  thickness.  Using  a  range  of  values  for  (3  similar  to  those 
chosen  for  a ,  the  lowest  error  (e  =  1.2)  is  obtained  when  a  =  0.3  and  [3  = 
0.7.  The  results  for  combined  flows  are  presented  in  Figure  6.  The  model 
compares  well  with  the  measured  combined  wave  and  current  friction 
factors,  and  the  average  wave  boundary  layer  thickness  of  6.2  cm 
determined  from  the  model  (a  =  0.3)  compares  well  with  the  measured 
value  of  6  cm. 

MM  also  conducted  experiments  for  pure  waves  (MM  1996a,  Table  1).  The 
roughness  elements  were  the  same  as  in  the  combined  flow  and  pure 
current  cases.  Given  that  MM  demonstrated  similar  roughness  for  waves  in 
the  presence  and  absence  of  currents  permits  an  additional  opportunity  to 
refine  the  closure  parameters  in  the  case  of  pure  waves.  Setting  the 
roughness  height  equal  to  the  average  obtained  by  MM  for  the  pure  current 
case,  and  using  their  experimental  input  wave  parameters,  friction  factor 
curves  are  generated  from  the  model  and  compared  to  their  data.  The 
results  are  shown  in  Figure  7.  The  lowest  error  (e  =  1.2)  for  the  friction 
factor  was  obtained  with  0=0.3  and  [3  =0.8.  The  average  wave  boundary 
layer  thickness  determined  from  the  model  was  6.0  cm.  When  (3  was  set 
equal  to  0.7,  the  lowest  error  for  the  friction  factor  still  occurred  with  a  = 
0.3,  but  the  modeled  wave  boundary  layer  thickness  had  a  mean  of  5.4  cm. 
In  both  the  combined  and  pure  wave  case,  a  consistent  result  emerges  in 
which  the  calibration  coefficients  maintain  similar  values.  It  must  be 
emphasized  that  the  suggested  values  for  the  closure  constants  are  only 
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valid  as  long  as  they  are  applied  to  the  stress  model  presented  above  (MM). 
Other  wave/current  bottom  boundary  layer  models  that  include  similar 
modifications  must  be  calibrated  before  they  are  used  in  applications. 


Figure  6.  Combined  wave  and  current  model  calibration  results  for  the 
closure  parameters  a . 
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Figure  7.  Pure  wave  model  calibration  results  for  the  closure  parameters  a  and  (3  . 


The  modifications  presented  above  are  not  without  a  theoretical  or 
empirical  basis.  The  first-order  correction  to  z1  leads  to  an  eddy  viscosity 
profile  that  is  similar  in  functional  form  to  expressions  developed  by 
Nielsen  (1992)  and  Sleath  (1991).  For  rough  oscillatory  flow  very  near  the 
bed,  Nielsen  (1992)  proposed  the  eddy  viscosity  K  =  0.004 Ab2kb2(>)  for 
Ab  /  kb  <  16 .  Similarly,  Sleath  (1991)  suggested  K  =  0.0025 Abkbco  in  the 
range  1<  A/K  <  120 .  Both  expressions  share  a  common  functional 
dependence,  namely  the  nonlinear  product  Abkbco ,  with  the  constraint  that 
c  +  d  =  2.  The  calibration  results  presented  above  are  formulated  in  terms 
of  an  eddy  viscosity  in  the  transition  layer  (z1  <  z  <  z2)  that  is  written  as 
K  =  ku,cwz1  ,  with  z1  =  alcw{l  +  1 6kb  /  Ab) .  Like  the  Nielsen  (1992)  and 
Sleath  (1991)  results,  this  modification  leads  to  an  eddy  viscosity  profile 
that  is  also  an  implicit  nonlinear  function  of  the  product  Abkb(o .  To 
illustrate,  a  pure  wave  in  which  the  roughness  is  large  enough  so  that 
z1  /  z0  <  1  is  considered.  In  this  case,  the  eddy  viscosity  becomes 
K  =  1 +  PK  /  A)  •  Noting  that  u*wm  /ub  =  Jfj2  along  with 

the  definition  of  ZU)m(=  kilwjti  /  co ) ,  the  eddy  viscosity  can  be  written  as 
K  =  CLK2fw  /  2 Abco(Ab  +  /3kb ) .  Expanding  this  expression,  for  example,  for 
constant  fw ,  gives  an  eddy  viscosity  that  is  proportional  to  Abkbco ,  which 
has  the  same  functional  form  as  Sleath  (1991).  Depending  on  the  definition 
of  fw ,  other  nonlinear  expressions  emerge.  As  an  example,  Kajiura  (1968) 
derived  a  friction  factor  of  the  form  fw  =  0.35(kb  /  Ab)2/3 .  Substitution  of 
this  expression  gives  K  =  0.175aK2(A^/3A:^/3Ci)  +  (3Ab3kb/3co) .  Although  this 
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functional  form  is  different  from  the  results  of  Nielsen  (1992)  or  Sleath 
(1991),  all  three  formulations  share  a  common  nonlinear  dependence  on 
the  same  external  parameters  kb,  Ab,  and  o) .  The  present  modification,  in 
that  the  explicit  dependence  on  kb  /  Ah  vanishes  in  the  limit  of  smoother 
bed  conditions,  is  somewhat  more  general  than  the  Nielsen  (1992)  and 
Sleath  (1991)  models  that  apply  only  for  kb  /  Ab  greater  than  approximately 
0.008. 

3.5  Speed  of  convergence  tests 

An  advantage  of  the  solution  algorithm  described  in  this  report  is  that  the 
nested  iteration  scheme  used  in  Styles  and  Glenn  (2000)  and  the  family  of 
Grant,  Madsen,  and  Glenn  models  (GM;  Glenn  and  Grant  1987)  can  be 
avoided.  Since  the  iterative  root-finding  algorithm  is  the  most 
computationally  expensive  step  in  the  solution  procedure,  it  is  of  great 
advantage  if  the  number  of  times  this  operation  must  be  executed  can  be 
substantially  reduced.  In  fact,  Keen  and  Glenn  (1994)  spent  considerable 
effort  to  optimize  the  initial  guess  for  the  friction  factor  and  other  variables 
to  speed  the  convergence  in  their  streamlined  version  of  the  GM  BBLM. 

To  illustrate  the  computational  advantage  of  the  present  approach  over  the 
Styles  and  Glenn  (2000)  model,  which  uses  the  same  eddy  viscosity  profile 
but  still  uses  a  nested  iteration  scheme,  results  of  a  speed  of  convergence 
test  are  presented.  Because  the  Styles  and  Glenn  (2000)  model  is  restricted 
to  a  much  narrower  range  of  wave,  current,  and  roughness  environments, 
the  input  parameters  represent  only  a  small  subset  of  the  full  capabilities  of 
the  stress  model  presented  here.  The  values  of  the  input  parameters, 
normalized  run-time,  and  total  number  of  iterations  are  listed  in  Table  2. 
Each  row  represents  10,000  independent  model  runs  with  identical  input 
and  initial  conditions.  The  run-time  was  recorded  for  each  run  and 
normalized  to  produce  the  numbers  listed  in  Table  2.  The  numbers  in 
parentheses  under  the  Styles  and  Glenn  (2000)  model  denote  the  maximum 
number  of  iterations  required  for  the  friction  factor  convergence.  The  other 
set  of  numbers  denote  the  number  of  iterations  required  for  u,c  convergence 
in  the  Styles  and  Glenn  (2000)  model  and  a  in  the  present  model.  The 
Styles  and  Glenn  (2000)  model  uses  the  secant  method  while  the  present 
model  uses  a  variation  of  Brent’s  root-finding  algorithm  (Atkinson  1989),  in 
which  the  iterations  are  performed  using  the  bisection  method  but 
convergence  is  checked  using  the  secant  method  after  only  a  few  iterations. 
A  tolerance  limit  of  icH,  or  a  0.01%  relative  error  between  the  previous  and 
present  iteration,  is  designated  to  established  convergence. 
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Table  2.  Speed  of  convergence  tests  comparing  the  present  method  with  the 
unstratified  version  of  the  Styles  and  Glenn  (2000)  wave  and  current  BBLM.  For  all 
model  runs,  (f>cw  =  0  and  Zr  =  100  cm.  The  first  three  rows  are  for  strong  waves  and 
currents  (SS),  the  middle  three  rows  are  for  strong  waves  and  weak  currents  (SW), 
and  the  last  three  rows  are  for  weak  waves  and  strong  currents  (WS).  The  last  two 
columns  list  normalized  run-time  (RT)  and  total  number  of  iterations  (N)  for  each 
method.  The  Styles  and  Glenn  (2000)  model  uses  a  nested  iteration  scheme.  The 
numbers  in  parentheses  indicate  the  maximum  number  of  iterations  for  the  inner 
loop,  which  usually  occurred  during  the  first  or  second  iteration  of  the  outer  loop. 


ut 

A 

“r 

K 

Present 

Styles  and  Glenn 

(cm/s) 

(cm) 

(cm/s) 

(cm) 

RT 

N 

RT 

N 

SSI 

50 

100 

20 

1.0 

1 

7 

3.6 

5(10) 

SS2 

50 

100 

20 

10 

1 

7 

2.8 

5(9) 

SS3 

50 

100 

20 

100 

1 

5 

2.9 

5(7) 

SW1 

50 

100 

1 

1.0 

1 

8 

4.4 

7(8) 

SW2 

50 

100 

1 

10 

1 

9 

2.6 

7(7) 

SW3 

50 

100 

1 

100 

1 

9 

2 

6(6) 

WS1 

1 

2 

50 

1.0 

1 

4 

21 

3(37) 

WS2 

1 

2 

50 

10 

1 

6 

7.3 

5(17) 

WS3 

1 

2 

50 

100 

1 

5 

7 

5(18) 

The  results  indicate  that  in  all  cases,  convergence  proceeds  with  fewer 
iterations,  is  at  least  twice  as  fast,  and  in  some  cases  an  order  of  magnitude 
faster,  thus  illustrating  the  greater  efficiency  of  the  bottom  stress 
algorithm  presented  here.  A  similar  speed  of  convergence  comparison  was 
performed  between  the  present  and  the  GM  model,  which  can  also  be 
formulated  without  a  nested  iteration  scheme  (Grant  and  Madsen  1986). 
The  results  were  similar  except  for  the  smoother  conditions  ( kb  =  1  or  10), 
in  which  case  the  present  formulation  usually  converged  20%  to  30% 
faster.  Since  the  present  model  uses  a  more  physically  realistic  eddy 
viscosity  profile,  has  a  correction  to  produce  accurate  estimates  of  the 
wave  boundary  layer  thickness  for  very  rough  beds,  and  is  more  efficient 
for  smoother  conditions  than  the  Grant  and  Madsen  (1986)  model,  it  is 
recommended  for  applications  in  which  estimates  of  the  near-bed  flow  and 
suspended  sediment  concentration  profiles  are  desired. 
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4  Summary 

The  authors  have  presented  a  fairly  robust  algorithm  to  compute  the 
enhanced  wave  and  current  boundary  shear  stress  components  for  an 
unstratified  bottom  boundary  layer.  The  stress  model  was  designed  for  a 
broad  range  of  wave  and  current  flows  and  was  formulated  without  the 
need  to  introduce  fictitious  currents  in  the  constant  stress  layer  to  obtain 
closure.  Instead,  the  governing  equations  for  the  wave  and  current  were 
reviewed  and  used  to  identify  important  velocity  and  length  scales  that 
could  characterize  the  flow  for  a  broad  range  of  wave  and  current 
conditions,  including  the  limiting  case  of  pure  waves  or  pure  currents. 
Systematic  non-dimensionalization  of  the  governing  equations  revealed 
three  important  internal  parameters:  e  =  U,c  /  U*cw,  /J.  =  U„wrn  /  u*cw  and 
O  =  Ub/  ILCW .  It  was  demonstrated  that  interpreting  the  functional 
dependence  of  e  ,  q ,  and  a  graphically  helped  to  illustrate  model  stability 
and  to  distinguish  the  effects  of  different  turbulence  closure  methods  (i.e., 
different  eddy  viscosity  formulations).  The  bed  shear  stress  was  most 
sensitive  to  the  form  of  the  eddy  viscosity  in  the  outer  wave  boundary  layer 
and  above  for  rough  flow  conditions. 

For  very  rough  conditions,  available  combined  wave  and  current  data  were 
utilized  to  refine  estimates  of  the  empirical  constant  a .  To  resolve  the 
discrepancy  between  past  formulations  that  produced  accurate  estimates 
of  the  friction  factor  but  underestimated  the  thickness  of  the  wave 
boundary  layer,  the  scale  heights  z1  and  8CW  were  modified  to  include  an 
explicit  dependence  on  the  relative  roughness.  This  introduced  an 
additional  closure  constant,  /3 ,  that  was  determined  experimentally  to  be 
approximately  0.7  for  combined  flows  and  0.8  for  pure  waves.  Further 
analysis  of  combined  wave  and  current  flows  over  very  rough  beds  in 
natural  flows  is  needed  before  a  definitive  value  can  be  prescribed  to 
model  the  constant  stress  portion  of  the  bottom  boundary  layer.  Until 
then,  it  is  suggested  that  a  =  0.3  and  p  =  0.7  for  applications  of  the  stress 
model  presented  here. 

Speed  of  convergence  tests  revealed  that  the  present  model  converged  in 
fewer  total  iterations  and  much  faster  than  the  Styles  and  Glenn  (2000) 
BBLM,  which  used  the  same  eddy  viscosity  profile.  Faster  convergence  was 
attributed  to  the  more  efficient  solution  method,  which  avoided  the  nested 
iteration  scheme  used  by  Styles  and  Glenn  (2000)  and  the  family  of 
Grant,  Madsen,  and  Glenn  models. 
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Appendix  A:  Derivation  of  Non-dimensional 
Wave  Shear  under  Strong  Waves 


Presented  here  is  a  derivation  of  rius  and  g  for  the  eddy  viscosity  in 
Equation  (2-5)  when  £0  is  greater  than  ^  or£2  .  Invoking  the  usual  linear 
and  boundary  layer  approximations,  Styles  and  Glenn  (2000)  present  the 
governing  equation  for  the  wave  within  the  wave  boundary  layer  for  the 
three-layer  eddy  viscosity  as 


iW 

iW 

iW 


d  r  dW  n 


dV 

d  *.  dw 


dV1  dt 

d  rdW 


=  0 


dV  d% 


0 


£2<£, 

£0<£<£l 


(A-5) 


where  W  =  uw  —  ub .  A  simple  harmonic  motion  (emt )  so  that  the  time 
dependence  is  separable  from  the  £  dependence  is  assumed.  With  the 
appropriate  boundary  and  matching  conditions  (Styles  and  Glenn  2000), 
the  solution  for  the  modulus  of  the  wave  is  written 


IG 


IG 


IG 


ub  +  G(Ker2^Je  +  iKei2^Je 

Ub  +  Ce^  +  De-”*  | 

ub  +  A{Ber2 ^  +  iBei2^)\ 


S2<Z, 

£l  <  %  <  ^2> 

£(><£<£ 


(A-6) 


where  A,  B,  C,  D,  and  G  are  complex  constants,  Ber,  Bei,  Ker,  and  Kei  are 
zero-order  Kelvin  functions,  and  m  =  ^Ji  /  .  For  convenience,  the  e""' 

term  is  dispensed  with  since  in  this  example  the  modulus  of  a  complex 
number  (f+  ig )  times  eicot  is  simply  the  modulus  of/+  ig.  The  values  of  the 
constants  can  be  found  in  Madsen  and  Wikramanayake  (1991)  and  reflect 
the  specific  condition  that  £0  <£i<£2  .  Based  on  the  governing  equation  of 
Equation  (A-5)  and  solution  of  Equation(A-6),  it  is  possible  to  extend  the 
wave  solution  to  include  cases  when  £i<£0<£2  and ^  < £2  < £0 . 
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Case  1:  ?1<?0<f2 


Using  the  eddy  viscosity  profile  given  in  Equation  (2-10),  the  governing 
equation  for  W  can  be  written  as 


— 0 

iw=i^=° 


Z2<Z, 

^0  ^  ^  ^2‘ 


Invoking  the  no-slip  condition  at  the  bed  and  given  that  the  solution 
smoothly  approaches  the  potential  flow  result  at  the  top  of  the  boundary 
layer,  the  corresponding  solution  for  the  wave  modulus  is 


uw  =  ub  +  G'(Ker2^Je  +  iKei2^Je  f2  <  f , 


uw  =  ub+C'emZ  +  D'e~mZ 


where 


C/_  ~ubL  ry  _  ~ubN 

PqL  +  MqN’  PqL  +  M0N 


n,_nC'P2  +  D'M2  ,  JC'P2-D'M2) 

V  'P  TS-fl) 


(A-10) 
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The  terms  in  Equations  (A-9)  and  (A-10)  are  defined  as  follows: 

M2  =  e~m^, 

=  e"^°,  P2  =  e171^, 

K2  =  Ker2^2/e  +  i  Kei2^2/e, 
d 


M0  —  e 
P0  =  er"b0 


KIP  =  ^(Ker2^J~e  +  i  Kei2^J~e) 

L  =  MJmK2  +  KfP),  N  =  PJmK2  +  KfP). 


?=?2 


(A-ll) 


Substituting  the  wave  solution  into  Equation  (2-7),  becomes 


‘  IPS 


m 


M0N  —  P0L 


MnN+PnL 


(A- 12) 


and 


I12  =  kPoTws. 


(A- 13) 


Case  2:  ?1<f2<?0 

For  this  case,  Equation  (2-13)  defines  the  eddy  viscosity  (Figure  1  [c])  so 
that  the  solution  for  the  wave  modulus  in  the  vicinity  of  £0  becomes 


1  Kev2^Jl  +  iKei2jf/e 
Ker2^jl  +  iKei2^J< 


(A- 14) 


Inserting  Equation  (A- 14)  into  Equation  (2-7)  gives 

d 


'  WS 


°^{Ker2jZJe  +  iKei2^J~e) 


?=?0 


Ker2^Je  +  iKei2^Je 


(A-15) 


and 


P 1  =  K^0eaTws 


(A- 16) 
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Appendix  B  -  Derivation  of  Non-dimensional 
Wave  Shear  under  Weak  Waves 


Presented  here  is  the  solution  for  and  p  while  neglecting  the  stress 
term  in  the  governing  equation  for  the  wave  above  £2 . 

Case  1:  Two-layer  eddy  viscosity  (z0  <  z1) 

For  the  two-layer  eddy  viscosity  presented  in  Figure  l  (d),  the  governing 
equation  for  W  is  similar  in  form  to  the  lower  two  layers  in  Equation  (A-5). 
Applying  the  appropriate  boundary  and  matching  conditions,  the  modulus 
of  the  wave  solution  becomes 


u 


IV 


U 


w 


ub  +  D"e  m?| 

ub  +  A\Ber2 ^  +  iBei2^) 


+B\Ker2^  +  iKei2^) 


^0  ^  ^  ^  ^i’ 


(B-l) 


where 


A'—  -U»N' 


B0N'  —  K0L'  ’ 


B 


/ _ 


~ubL 


K0L'  —  B0N' 


(B-2) 


M1(l  — m) 


(B-3) 


and 
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The  terms  in  Equations  (B-2)  and  (B-3)  are  defined  as  follows: 

K0  =  Ker2^f0  +  i  Kei2^,  B0  =  Ber2^f0  +  i  Bei2^, 

K±  =  Kev2 ^ + i  Kei2^,  B1  =  Bev2 ^ + i  Bei2^, 

d 


K?  =  ^(jfei-2^ + i  Kei2, ^ 


Bi»  =  -^[Ber2^f1+iBei2^f1 


?=?1 


£=fi 


L'  =  mB1  +  B{1] ,  N'  =  mK1  +  K{1] , 

M'  =  e~mk. 

Substituting  the  modulus  into  Equation  (2-7)  yields 


N'BP 

0 

vl 

J 

BqN'  —  KqL' 

K0L'-B0N'\ 

and 

li2  =  K%0(lTws. 

Case  2:  Two-layer  eddy  viscosity  (z0  >  z1 ) 

If  £0  is  greater  than  ^ ,  the  wave  modulus  becomes 


U„ 


Uu 


£1  ^  ^0  ^ 

With  the  aid  of  Equation  (2-7),  takes  on  the  simple  form 


r ws  =  m 


and 


li2  =  K^oTws  =  K^O. 


(B-4) 


(B-5) 


(B-6) 


(B-7) 


(B-8) 


(B-9) 


7.  (concluded) 
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